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Abstract 

We consider a model for the polymerization (fragmentation) process involved in infectious 
prion self-replication and study both its dynamics and non-zero steady state. We address 
several issues. Firstly, we give conditions leading to size repartitions of PrPsc aggregates 
that exhibit bimodal distributions, as indicated by recent experimental studies of prion ag- 
gregates distribution [19]. This is achieved by a choice of coefficients in the model that are 
not constant, thus extending a previous study of the nucleated polymerization model [16] . 
Surprisingly, conditions for bimodality do not seem to be the most favourable for prion repli- 
cation. Secondly, we show stability results for this steady state for general coefficients where 
reduction to a system of differential equations is not possible. We use a duality method based 
on recent ideas developed for population models. These results underline the potential influ- 
ence of the amyloid precursor production rate in promoting amyloidogcnic diseases. Finally, 
we numerically investigate the influence of different parameters of the model on PrPsc accu- 
mulation kinetics, in the aim to study specific features of prion strains. This study suggests 
that PrPsc aggregate size distribution could be a signature of a strain in a given host and 
a constraint during the adaption mechanism of the species barrier overcoming, that open 
experimental perspectives for prion strain investigation. 



1 Introduction 

Transmissible spongiform encephalopathies (TSE) are fatal, infectious, neurodegenerative dis- 
eases. They include bovine spongiform encephalopathies (BSE) in cattle, scrapie in sheep and 
Creutzfeldt- Jakob disease (CJD) in human pp. The infectious agents responsible for disease 
transmission, known as prions, present some unusual biological properties (as a high resistance 
to inactivation by heat or radiation). According to the "protein-only hypothesis", prions may 
consist in a misfolded protein (called PrPsc), without any nucleic acid. This hypothesis suggests 
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that PrPsc replicates in a self-propagating process, by converting the normal form of PrP (called 
PrPc) into PrPsc [2]. Many evidences are in favour of an autocatalytic replication of PrPsc, as 
the generation of infectivity from recombinant proteins [1] or the use of in vitro PrPsc conversion 
systems, such as the protein misfolding cyclic amplification (PMCA) technique [5]. 

However, the precise mechanism of conversion remains unclear. Moreover, prion infectious 
agent can exist under different strains, characterized by their incubation period and their lesion 
profile in brains [3]. In the framework of the protein-only hypothesis, it is supposed that strain 
diversity is supported by various conformation states of PrPsc, that leads to various biological 
and biochemical properties [6j[7]. A critical challenge of prion biology consists in elucidating the 
mechanism of conversion of PrPc into PrPsc, and therefore how a diversity of strains may exist 
in the same host (expressing the same PrP molecule). 

To investigate the conversion of PrPc in PrPsc, many relevant mathematical modeling of prion 
replication have been proposed [H [91 Uni E] • Major aim of these modelling is to demonstrate that 
essential features of prion disease can be explained by purely physico-chemical mechanisms, as 
supposed by the protein-only hypothesis. In addition, mathematical modeling allows to study the 
effect of every elementary process in a separate manner [12] , what is difficult to do experimentally. 

The early proposed model is the heterodimer one. It is based on the conformational change of 
PrPc into PrPsc after the formation of a heterodimeric complex (PrPc + PrPsc — > PrPc*PrPsc 
— > PrPsc*PrPsc — > 2 PrPsc). This model does not take into account the aggregation of PrPsc, 
and thereby fails to explain the association between infectivity and aggregated PrP. Some other 
mechanisms have been proposed, which are interested in PrP aggregation [9| 1 10| . IT3l IT^l I15j . Based 
on fibrilar aggregation, the model which seems by now broadly accepted is the one of nucleated 
polymerization (9J 116] . In this approach, PrPsc is considered to be a polymeric form of PrPc. 
Polymers can lengthen by addition of PrPc monomers, and they can replicate by splitting into 
smaller fragments. It is worth noting that deterministic [U [TBI [T7] as well as stochastic [121 E] 
simulations of this model lead to an unimodal size distribution of PrP aggregates, which seems 
to be quite insensitive to small variations of parameters [12]. Greer et al. [18] recently improved 
the model and include a mean saturation effect by the whole population of polymers onto the 
lenghtening process (called general incidence), and polymer joining (through a Smoluchowski 
coagulation equation). In these models, each aggregate has the same behaviour, regardless to its 
size. However, recent experimental analysis of relation between infectivity and size distribution 
of PrPsc aggregates (for PrPsc purified from infected brain [19] or for PrPsc produced by PMCA 
[201121] ) contradicts this uniform behaviour of PrPsc aggregates. In addition, Silveira et al. shows 
that bimodal distribution of the polymer size is more likely to occur within the real process [19] . 

The goal of this present study is to better understand the unexpected experimental size 
distribution of prion aggregates in brain, and to investigate the potential implication of this size 
distribution in the strain phenomenon. To do so, we generalize Masel and coauthors' model so 
as to take into account aggregate size-dependent parameters. 

The paper is organized as follows: in section [2] we recall and review the model of Masel et al. 
and its continuous version [16] which is going to be used. We introduce our main improvement, 
namely a size-dependent lenghtening factor related to nonuniform infectivity rate; and we discuss 
the eigenvalue problem which is a key tool to analyze the model. In section [3] we investigate 
heuristically and numerically the formation of a bimodal distribution of polymers. In section 0] 
we prove the stability of the zero steady state in the disease free regime (generalizing partially 
a result in [16] ) . Finally, in section [5] we study numerically the influence of different parameters 
on the dynamics of our model, in a prion strain perspective. 
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2 The continuum model 



The following set of coupled differential equations has been introduced by Masel et al. [9] in 
order to model the polymerization (aggregation, fragmentation) process involved in infectious 
prion self-replication. It describes the dynamics of the quantity of PrPc, V(t), coupled with the 
evolution of aggregates of PrPsc which have size i, Ui(t), 



dV_ 
~dt 

du 



71()-1 

A - ~/V - tVU + 2/3 X/ M J' 

i=l j>i 
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The index no denotes the minimal size of PrPsc polymers. The quantity U(t) = ^2ui(t) is 
the total amount of prion aggregates. The constant parameters A, 7, r, (3, u are respectively the 
basal synthesis rate of PrPc, the degradation rate of PrPc, the conversion rate of PrPc into 
PrPsc (autocatalytic process following the mass action law), the fragmentation coefficient and 
the degradation rate of PrPsc. 

Analysis is simpler in the framework of continuous size of prions because analytical tools can 
serve to find simpler formulations. Accordingly, Greer et al. introduce a continuous version of 
([1]) where the variable x £ (0, +00) denotes the size of aggregates and replaces the index i £ N. 
Asymptotic derivations of continuous models from discrete models can be found in [22^ 123]. The 
continuous model reads, with possibly nonconstant coefficients, 



dV(t) 
dt 
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t(x)u(x, t) dx 



■t'o 
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x / P{y)K(x,y)u(y,t)dy dx, 



u(x, t) = -V(t) — (r(x)u(x, t)) - \n(x) + P(x)]u(x, t) + 2 / P(y)it(x, y) u(y, t) dy, 



u(x ,t) = 0, 



together with appropriate initial conditions. This is a well established family of models used 
for describing aggregation, fragmentation and possibly coagulation in polymers and also size 
structured cell dynamics with finite resources [Ml ESI [261 12Z]- Well-posedness, in the class of 
weak solutions, can be found in [28, 129] . 

d 

The transport term V(t) — (r(x)u(x,t)) accounts for the growth in size of polymers: their 

size grows with the speed V(t)r(x), proportional to the available PrPc molecules V(t), with an 
aggregation ability depending on the size of the polymer (a conceivable hypothesis being that 
their size confers them a peculiar geometry affecting the autocatalytic process). The fragmenta- 
tion rate, for a polymer of size y, is (5{y) > 0. The repartition of the two fragments of (smaller) 
sizes x and y — x is given by k(x, y) > 0. It should thus satisfy the two usual laws [30] expressing 
that the number of fragments increases but with constant total molecular mass (recall the factor 
2 in the right hand side of ([2])) 



k(x, y)dx = 1 , 



x k(x, y)dx = — 



(3) 



As usual for size-structured population models [25, 30J, this model may incorporate a minimal 
size of infectious PrPres aggregates xq > 0, whose value remains unknown. Experimentally, no 
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monomer of PrPsc has been isolated yet. In addition, small aggregates have been shown not to 
be infectious [19] (even-though they have to enter the actual modeling) and thus the assumption 
of a critical size of nucleation xo has been used. However, the continuous model holds under 
the assumption that the monomer size can be neglected as opposed to possibly numerous large 
polymers (such that xq ~ 0). Therefore we keep xq in the following but we sometimes neglect it 
in order to simplify the presentation. 

Assumptions on the coefficients. In its entire generality, model ([2]) is rather difficult to 
attack even though some qualitative behaviors can be described as in section [H We aim to 
reduce the complexity of this system in order to extract some relevant information as the bimodal 
distribution motivated by the recent work of Siveira et al. [19] . A simple choice for the coefficients 
is as follows [16J: the fragmentation rate is taken proportional to the fragment size and the 
degradation rate does not depend of the size, that is 

0{y) = 00 V, fJ>(x) = Ho, (4) 

and the probability distribution of fragments of size x is chosen to be uniform with respect to 
the length of the splitted polymers of size y: 

, , |0 if y < x or y < x ... 
k{x, y) = < (5) 
I \ jy if y > xo and < x < y . 

Last, but not least, the transition rate t(x) may also be assumed to be constant: 

t(x) = Tq. (6) 

In this situation the possible time asymptotic behaviors have been completely classified [91 [16] 
as well as the stable distribution of aggregates (see Figure E^left) and the paragraph "Constant 
coefficients" at the end of this section). When referring to the model of Greer et al. with a 
constant rate To, we will denote it by the 'constant coefficients' case. 

With this respect, the main purpose of our work is to take into account different shapes for 
the function t(x) in order to fit the recent experimental data. Also we perform a stability anal- 
ysis which does not assume particular coefficients as in (HH6|) . At this stage of knowledge of the 
biochemical process, the available microscopic experimental data are insufficient to investigate 
different fragmentation laws in details. Indeed, we are not able to dissociate the three elementary 
processes (degradation/sequestration, splitting and polymerization) implicated in prion replica- 
tion. We have made the choice of varying the rate of infection t(x) as a first step, whereas 
Silveira et al. clearly indicate that large polymers are more stable than small ones (this could 
be assumed as a reduced degradation rate or a reduced fragmentation factor for large sizes). 

The system ([2]) keeps an important biochemical property, that is the prion molecules are 
properly transfered from one configuration to another (inducing no loss of mass during frag- 
mentation or polymerization). This enhances the following macroscopic law involving the total 
quantity of polymers U and the mean length of aggregates P: 

roo poo 

U(t) = / u(x,t) dx , P(t) = / xu(x,t) dx . 

In fact equation ([2]) yields, 

d f'°° 

— (V(t) + P(t)) = A - -yV(t) - / ii(x)xu(x,t)dx . (7) 
at J xo 
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Figure 1: Experimental analysis of PrPsc aggregates. (Data have been kindly provided by Silveira. 
For details, see [19] ). Size-distribution of PrPsc aggregates in a whole infected hamster brain 
(Solid line) and relative specific infectivity of each fraction containing PrPsc aggregates with 
respect to their size (dotted line). The former should correspond to the quantity xu(x,t) within 
our formalism, and the latter represents a bell-shaped converting rate being t(x) in model ([2]). 
Abscissa = fraction number; ordinate = PrP (in ng) for the solid line, and relative specific 
infectivity (in arbitrary units) for the dotted line. Note that the scale of the polymer size x in 
this graph is not likely to be linear. However we can hardly evaluate the precise scale because 
of the limitation of measurement experiments. 
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Stationary polymer distribution. Contrary to the one-peak distribution of infectious prion 
proteins depicted in Figure [3jleft) under the 'constant coefficients' assumptions (jlHS]), experi- 
ments indicate that this distribution is rather bimodal (Figure [1]). It comes out that bimodal 
configurations could be the result of a non-constant converting factor t(x) (having a bell shape 
for instance). As we are interested mainly in long-time dynamics resulting in a balance between 
polymerization and fragmentation, we will often stress out the so-called 'Stable Size Distribu- 
tions' (see [3D]) that are the stable steady states of system ([2]), i.e. satisfying 



7 + / t(x)uoo(x) dx 



rxo 

A + 2y x\ (3(y)K(x,y)u 00 (y)dydx, 



d 

^oo^(t(x)uoo(x)) + \fi(x) + P{x)]uoo(x) = 2 j (3{y) k(x , y) u^iy) dy , 



(8) 



with sufficient decay conditions at infinity. There are in general two steady states, the first one 
corresponds to 'healthy cells', and we refer to it as the 'zero steady state': 

V = -, u = 0. (9) 

7 

A second possible steady state, denoted by (Voodoo), corresponds to the infection regime, and can 
be understood as follows. For a given V, the fragmentation equation has a dominant eigenvalue 
A(V) (growth rate), associated to a nonnegative eigenvector. In other words there is a unique 
solution to 

V— (t(x)U(V; x)) + \fi(x) + 0{x)]U{V- x) - 2 J /%)«(*, y)U(V; y) dy = A(V) U{V; x), 

U(V;x ) = 0, U(V;x)>0, f U(V; x) dx = 1 . 

(10) 

Such solutions have been shown to exist, and the integrability condition implies a much faster 
decay at infinity, see [251 HSU EI] . 

Remark 1 . Intermediate asymptotics of ([5]) and particularly early stages of the infection ( expo- 
nential phase) are also considered in the sequel (Section^. Interestingly, the distribution shape 
corresponding to the exponential phase (Figure\B(left, full line)) matches quantitatively with the 
experimental data of Figured Indeed we observe that large polymers are prevalent, corresponding 
to a level of PrPc close to the healthy state V V (see also Section^. 

The equilibrium distribution Uoo arises as an eigenvector associated to the eigenvalue A(V^o) = 
0. In fact, this characterizes both the ground level of PrPc, V^,, and the shape of the polymer 
distribution ^(V^; z), up to a constant factor (the total number of polymers). This missing 
factor is determined thanks to the first equation of ([8]) as we show it later. Interestingly enough, 
the value does not depend on the differential equation driving V(t). In particular it does not 
depend on A and 7. 

The decay condition at infinity ensures that the following integrations by parts can be justified 
(we assume xq = below for the sake of simplicity). Integrating (llOj) successively against 1 and 
x, and using ([3|), gives respectively 

/p pCsO pCsO p 

n(x)U(x) dx + J (3{x)U{x) dx - 2 J J f3(y)n(x,y)U(y) dy = A J U{x)dx , 
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Figure 2: Eigenvalue problem for the microscopic distribution, (left) Several eigen- 
functions are plotted, for V respectively above (full line), below (dotted line), and close to 
(dashed line) the equilibrium value VJx>. This is the 'constant coefficient' configuration except for 
the bell-like converting rate t(x). Coefficients' values are given in Section [5j (right) Numerical 
computation of the eigenvalue function A(V) for a bell-like function t(x) (dashed line) compared 
to the 'constant coefficients case', where A(V) = (jlq — y/ro^oV (full line). 

and 

/r /'OO roo r 

t(x)U(x) dx + / (xfi(x) + x(3{x))U{x) dx — 2 / x n(x, y)[3(y)U(y) dy = A / xU{x)dx 

As a direct consequence we obtain 

a f\r\ [r,, M s\,j(v \ A ~Xl ^ x)dx + f x^{x)U{V- x) dx 

A(V) = J frix) - 0(x)]U(V; x) dx J xU (V;x)dx " (U) 

Macroscopic (non zero) steady state. Back to ([8]), the steady state of interest, if it exists, 
is given according to the rule 

A(F oo )=0, u 0O (x) = g oo U(y 0O ]x) , (12) 

where g denotes the total number of polymers. There do not always exist non-trivial steady 
states. To see it, consider again the case xq = 0. Then, according to the first equation of system 
(|SJ), we have 



or, equivalently 



1 + Qoo j t(x)U(Voo]x) dx 



A, 



> . 



/ T(x)U{V 00 ;x)dx 

This points out a constraint for the non trivial steady state (V^-Uqc) to exist, namely it is 
required that 

7 Foo<A, (or VooKV). (13) 

We tackle in this issue the nonlinear (in)stability of the trivial steady state (V, 0) in section 
d] for general coefficients. In particular, we show that stability of the trivial steady state when 
V < Voo follows directly from the property that 

A(V) is a decreasing function. (14) 



Constant coefficients. The case of 'constant coefficients' (J^HSJ) has been completely under- 
stood by Greer et al. [16J. We have the opportunity to recall their results under the viewpoint 
of the eigenvalue A in (fTUj) . This eigenvalue can be explicitly computed from (fTT|) which yields 

A(V0- w = -ft/,M(V,, )<fe ji -y i)|fa . 
Eliminating the quantity J x U(V; x) dx = %p we obtain 

A(V0 = Mo - V^A^ 7 • (15) 

Notice that it is a decreasing function of V (this property is in fact crucial in the subsequent 
analysis). We also point out for further use the following consequence for the non-zero steady 
state: 

J x u oc (x)dx _ /io ^ 
/ Uoo(x) dx fio 

Following the constraint (|13p . a non trivial steady state exists if and only if 7 < Xtq. Still from 
|16j . we know that it is globally asymtotically stable. On the contrary, when it does not exist, 
the zero steady state is globally asymptotically stable. 

3 Bimodal fragment distribution of PrP 

Recent experimental data obtained in [19J suggest that PrPsc aggregates have different infectious 
abilities, according to their size (see Figure [1]). Although infectious ability actually seems to 
be a competition between the three processes involved in the model (namely fragmentation, 
polymerization and degradation) one simple way to improve the model compared to [16J is to 
introduce a single variable coefficient. Results from [19] show different in vitro abilities to convert 
PrPc into PrPsc, arguing in favour of a non-constant polymerization rate. 

Therefore, in this section we consider the system ([8]) with the choice fj,(x) = /io, f3{x) = (3qx 
and n given by ©, and a non-constant polymerization rate r = t(x). We study, by means of 
numerical experiments and analytical methods, how bimodal fragment distributions may occur. 

Direct numerical simulations. Our first observation is that a reasonable bell shaped con- 
verting factor t{x) induces a bimodal stationary distribution (see Figure [2] for instance). This 
is supported by various simulations of the model dH). This numerical result leads to investigate 
analytically which conditions are needed for a bimodal distribution to appear. This is our next 
step. 

A differential equation formulation. As a first step, we present a simpler differential for- 
mulation of the stationary problem ([8]) , where we strongly use our particular choice of coefficients 
(all are constants except t(x)). 

Differentiating equation ([8]), we find the following second order differential equation 

V 00 (t(x)u 00 (x)) xx + (""00(2) (no + 0qx) ) x + 2f3 u oo (x) = . (17) 

This differential equation is complemented by the boundary conditions 

/•oo 

Uoo(x ) = 0, Vooiruoo)^^ = 2(3 / Uoo(y) • 
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Figure 3: The splitting effect, (left) A typical stationary distribution for model ([2]) with 
'constant coefficients' according to (HHS]). This unimodal shape is given explicitly by a suitable 
dilation of the profile $(r) = (r + r 2 /2)exp(— r — r 2 /2). (right) For a suitable bell-shaped 
converting rate r(x) (dotted line), the non-zero steady state solution to ([8]) presents a bimodal 
shape (solid line). As compared to the constant converting rate case t(x) = tq (left), the peak 
of r acts so as to split the distribution into a bimodal one. Abscissa = PrPsc aggregates size; 
Ordinate = t(x) either u^x) in arbitrary units 

Recall also the additional relation derived from ([8]), after testing against x, 

Voo / t(x)u oc (x) dx = fi / xu oc (x)dx. 

Jxo Jxo 

The equation (|17p is reminiscent of the standard one-dimensional Fokker-Planck equation. 
In fact, if we drop the last term of order 0, it becomes 

V 00 (t(x)(u 00 (x)) x ) x + (tp x {x)u 00 (x)) x = , 

that is, a inhomogeneous diffusion with an effective potential given by ip{x) = V 00 t(x) + fiox + 
T^Pox 2 . Hence the contribution of a non constant t(x) might rise a double well potential on i/j(x), 
and thus be responsible for the formation of a two-peaked solution. This holds, by perturbation 
analysis, at least for (3q small. But we can give a more accurate estimate. 

Condition for a bimodal equilibrium distribution. We now give a necessary condition 
for a bimodal distribution. We evaluate the derivatives in the microscopic equilibrium equation 

dm), 

V 00 {t(x){u 00 ) xx {x) + 2t x (x){u 00 ) x {x) + t xx (x)u 00 (x)) = (-Uoo (no + Pox)) x - 2f3 u oo (x). 

A bimodal length distribution has a critical point where it is convex and therefore we are 
looking for a point x* such that {uoo) x {x*) = and (uoo) xx (x*) > 0. A necessary condition on 
t{x) for the existence of such a point is therefore 

inf Voo t xx (x) < -S/3q. 

x>0 

The above condition is indeed intricated because Voo itself depends on t(x). But it gives a first 
insight of the desirable conditions, meaning that a peak which is concave enough is required. As 
we will see it below, this condition is not sufficient. The locus of this strongly concave 'peak' is 
also important 
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Figure 4: Influence of the peak LOCUS. We superpose successively the function t[x) (full 
lines) and the corresponding equilibrium distribution (crossed lines) for translated bells. In each 
of these pictures, the center of mass of the polymer distribution is 1. We can visualize the 
splitting effect of the potential t{x) very well. 



Locus of the peak of t(x): heuristics. For coefficients like n(x) = fiQ and (3{x) = (3qx, we 
know from (|ll[h that the center of mass at the equilibrium configuration has the value /xq/A)- 
Departing with To constant, the unimodal distribution is split into a bimodal one under the 
action of t{x). The further is the peak of the bell-like function t{x) from /io//3o> the weaker is 
the splitting effect on the unimodal distribution. This is in accordance with numerical results 
(see Figure H]) . 

4 Stability 

We have seen previously that there are two possible steady states, setting the alternative between 
a disease free and an infected system and depending upon the prion production and degradation 
rates A, 7. For 'constant coefficients', Greer et al. [16J could study their stability using that the 
system can be reduced to a differential system on (V(t), U(t), P(t)) (see Section [2|). Here, we 
investigate the same question for general coefficients. Our main assumption is some monotonicity 
of the eigenvalue in problem (jlOp . which we have checked for 'constant coefficients', see formula 
(fT5j) . We choose xq = throughout this section. 

4.1 Stability of the zero state for V < 

We first tackle the stability of the disease-free steady state. Recall the notation V = — for the 
'zero steady state' from Q. 

Theorem 1 (Local stability). Assume that A(-) defined in ( fiflj) is decreasing and that V < Voo, 
where is defined by A(V r oc ) = according to HM) . Then, in equation (TJJ), the zero steady state 
V = V , u = is locally nonlinearly stable. 

We recall that in the case at hand (V < Voo), there does not exist a non-zero steady state 
because (fT3|) cannot be fulfilled with g > 0. We give in the subsection 14.31 below a class of 
coefficients for which we can compute the eigenvalue A(-) and it is indeed decreasing. 
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Proof. The eigenvalue A(-) is assumed to be a decreasing function of V. So, the condition V < 
ensures that A(V) > 0. We consider a perturbation of the ground state V(t) = V + V(t) and 
u(x,t) = + u(x,t) (note that u is nonnegative, whereas V has no sign). 
The nonlinear system for this perturbation writes 



d_ 

at 



dV(t) 
dt 

u(x, t) 



-V 



7 + / t(x)u(x, t) dx 



V / t(x)u(x, t) dx, 
Jo 

V^(T(x)u(x,t)) -V(t)^-(r(x)u(x,t)) -l^(x)+P(x)]u(x,t) (18) 

/•oo 

+2 / P(y)n(x,y)u(y,t)dy. 



Following the duality method in [32], we introduce the adjoint eigenvector for the eigenvalue 
V, namely Tp > given by the equation 

- Vr{x)^-lp{x) + [fM(x) + 0(x)]lp(x) = 2 I 0(x)K(y, x)<p{y) dy + A(F)^(x) . (19) 



We now introduce constants Ki, K2 such that 



< Knp(x), t(x) < K 2 (p(x) . 



(20) 



This is possible because Tp grows linearly at infinity according to general abstract properties 
proved in [251 113 EI] • See again subsection 14.31 for examples. 

We test the adjoint equation against the equation on u in (fT8l) and obtain 



d 



— I u(x, t)(p(x) dx = — A(V) / u(x, t)(f(x) dx + V(t) j ( — f>{x) ) t(x)u(x, t) dx . 



d 



dt 



On the other hand, multiplying the first differential equation in (|18p by the sign of V, we get: 



><-iki 



7 + / t(x)u(x, t) dx 



+ V t(x)u(x, t) dx . 



We obtain, choosing a large enough such that SA(V) 



d 



K 2 V 



>0, 



— [a j uZp + \V\) < -ak(V) j ulp + aKi\V{t)\ j Tpu-^\V\-\V\ j tu + K 2 V j uTp 



< - min A(V) 



K 2 V 



a 



7 [a / u7p + \ V\ + aKi\V(t)\ / Tpu 



A", 



< — min(5, 7) ( a j u<p + \ V\ ) H — — ( a j uip + \ V 



From this differential equation we conclude that, when ya f utp + \V\) is initially small enough, 
then the right hand side is negative. Therefore it decays for all times with the asymptotic 
exponential rate 

(a J uTB+\V^j <Ce- min ^ . 



□ 
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We can also state the following variant of Theorem [T] 

Theorem 2 (Global stability). Additionally to the hypotheses of TheoremUl assume thatr(x) > 
kTp(x) for some constant k > and V / A(V) is small enough compared to kj{K\K2)t. Then, in 
equation |2]), the zero steady state V = —, u = is globally nonlinearly stable. 

By opposition to those of Theorem [TJ the assumptions of this Theorem are difficult to check 
directly because the coefficients are intricated here. They mean that we are close to the 'constant 
coefficients case' because Tp = 1 in this case and we can choose K\ = 0. 

Proof. With these additional assumptions, we may keep one negative term in the last computa- 
tion and arrive to 



d_ 

It 



(a j wp + lV^ < -aA(V) J uTp + aKx\V(t)\ j Tpu - j\V\ - \V\ j tu + K 2 V j uTp 

< -min (a(V) - ^-.7) ( a J u^+\V\^j +{aK 1 -k)\V(t)\ J Tpu 

< -min ^A(F) - ^p7^ (ajwp + \V\J . 

This last inequality will hold if we can find a > j^m such that aK\ < k, which is precisely 
our smallness assumption. Then we have exponential decay of the solution. □ 

4.2 Unstability of the zero state for V > 

The same kind of method allows for another type of results initiated in [32], namely that solutions 
cannot go extinct when another steady state exists. Again we choose xo = in this subsection. 
We have, recalling that Tp is defined in (|19p . the 

Theorem 3 (Unstability). Assume that A(-) defined in U0\) is decreasing, that for some constant 
K > 0, t(x) d x Tp(x) < KTp(x), that V(0) < V and that < V. Then, in equation |lj), the zero 
steady state V = V, u = is locally unstable. 

Proof. Following the proof of Theorem [TJ we have on the one hand 

d — — f°° 

— (V- V{t)) + 7 (V - V(t)) = V{t) / t(x)u(x, t) dx > 0. 

dt Jq 

Therefore we have V(t) < V for all t > 0. 

On the other hand, we can again combine the equations (|19p for Tp and the equation on u in 
([2]), to obtain 

^ roo roo Q roc 

— / u(x, t)Tp(x) dx = (V(t) — V) j r(x)—Tp(x) u(x,t) dx — A(V) / u(x,t)Tp(x) dx. 

dt Jq Jq OX Jq 

We define the quantity (sum of two positive terms) 

("OO 

Q(t) = / u(x,t)Tp(x)dx + V -V(t). 
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As a consequence of the above calculations and of the hypothesis r d x (p < Kip, we have therefore, 
d 



d -Q(t) > -K(V-V(t))l Tp(x) u(x,t)dx- A(V) I u(x,t)Tp(x)dx 



Jo 

oo 



> ki[l-k 2 (V -V(t))] / u(x,t)lp(x)dx, 

Jo 

for some constants k{ > 0. Indeed, in this case when VJx, < V, we know that A(V) < 0. 

We can conclude from this inequality that is unstable. Indeed, whenever Q(t) becomes 
small enough, then V — V becomes also small enough so that Q(t) increases and the solution 
gets away from 0. □ 

4.3 A class of examples 

In order to clarify the assumptions and properties stated before, we give a class of coefficients in 
equation ([2]) where the different eigenelements can be explicitly computed. In each case we will 
see that A(V) s a decreasing function of V, and that (|20p reduces to the fact that r(x) is bounded. 

We first recall from [25\ 131] . the case r = To, fi = /xq an d (3 = (3\ constant. We obtain the 
solution to (f2~T1) below, 

A.(V) = Pi, (fix) = constant. 

The eigenlements are usually difficult to evaluate in the direct problem (|10p but can be easily 
computed in the adjoint equation, recall (fT9l) . 

-Vr{x)^tp{x) + \n{x)+p{x)](p{x) =2^ 0(x) K (y,x)ip(y) dy + A(V)<p(x) . (21) 

Indeed, searching for an affine solution (p(x) = 1 + x/L, we find, using the structure properties 
of k(x, y) in ([3]), 

= *^ +£(*) + (A - M (z))(l + ! 

Hence a more general class where one can compute A(^) is when (3(x) = (3q x, and tq, are 
constant. We arrive at 

-f3 L 2 + r V = 0. 



Therefore we obtain L = y/ro/fio and A(V) = fio — f3oL(V), as in Section [2l Observe that A(V) 
is a decreasing function of V as used in Theorem [H 

We leave to the reader to check that A(V) is also decreasing when j3(x) = j3\ + fio x, r(x) = 
To + T\ x and is constant. Indeed solving the two equations for (p(x) = 1 + x/L, we find 

A- fi + — +/3 1 = 0, — h -j- + fa = 0. 

And thus Z = A — /xo solves 

(Z + /3 1 )(Z + yr 1 ) = yr /3o, ^ < -Ft x . 
This last constraint comes from the condition L > 0. 

More generally, notice that it always holds true that 

A(0) = /i , A(V) < fi . (22) 
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To prove it, we just multiply (fTOj) by x and integrate. Notice that the problem is then degenerate 
and the eigenfunction is singular at x = 0. 

On the other hand, still for the same reason, if min > 0, we have 



and then there is at least one non-zero steady state. 

5 Temporal dynamics of the aggregates size distribution in a 
biological context 

As we highlighted previously, recent experimental work by Weber et al [20|, [2T] and by Silveira et 
al [19] suggests that PrPsc aggregates are differentially infectious with respect to their size. The 
former also shows that the PrPsc size distribution in a whole infected brain is bimodal rather than 
unimodal. Based on the nucleated polymerization model for prion growth, we have shown that 
bimodal stationary distribution can occur for a size-dependent conversion rate. However, during 
the time course of incubation period and clinical stage of prion diseases, PrPsc accumulation 
seems to follow an exponential dynamic until the death. Apart from hemizygote mice (PrP+ 
/0), no steady state of PrPsc accumulation was observed in the brain during all the disease 
\33\ [34] . Consequently, we focus in this section on the dynamic of PrPsc accumulation in the 
early stages of the infection named 'exponential phase'. Our approach is based both on the 
eigenvalue problem and on direct numerical simulations of the temporal dynamics. 

Numerics. Parameters of the nucleated polymerization model for prion growth have been 
estimated for the 'constant parameters model' from experimental data and are available in the 
literature [HI [161 H2] ■ The parameter values used in the sequel have been quoted from Rubenstein 
et al. [12]. Unless explicitly mentioned, they are: A = 2400 per day, 7 = 4 per day, fiQ = .05 per 
day, and = 0.03 per day. The conversion function t[x) is the sum of a basal rate tq = .001 
and a gaussian bell centered on to: 



with a magnitude A to be chosen by several orders of magnitude above tq (from biological 
evidence on different specific converting activities in [19]). The simulations assume an initial 
PrPc population V(0) = V = A/7 (corresponding to the healthy steady state) and an initial 
PrPsc population given by u(x, 0) = 0.5x 2 /(l + x 4 ), which is a small perturbation of the zero 
steady state indeed. 

Dynamics of PrPsc accumulation and Eigenvalue problem. The eigenvalue problem 
can be used to measure the dynamics of PrPsc aggregates accumulation during the exponential 
growth. We assume below that V > Voo, that A(V) is a decreasing function and that instability 
of the healthy state holds true (see Section HJ) . The practical choice of coefficients is mentioned 
above. 

Assuming that V(t) remains close enough to V (this is exponential growth phase in a lin- 
ear regime), then the polymerization behaves following a linear problem. Actually the second 
equation in (|2J) is decoupled from the first one at the first order of approximation. The domi- 
nant eigenvalue —A(V) thus measures the exponential growth of the PrPsc total population (see 
|25j for a mathematical formulation of this fact using the generalized relative entropy). As a 



A(oo) 



= —00 



(23) 




(24) 
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by-product, we also learn that the renormalized distribution tends to align along the eigenvector 
U(V,x). One noticeable difference with the stationary case studied in Section [2] is that the cen- 
ter of mass of U(V,x) is translated to the right as opposed to lify^x) (for which it is Ho/Po), 
according to (fTT|) : 

0<-^= f (x-^-)u(V,x)dx. 



A, J V /V 

Observe that for the parameters picked up in the literature, the difference is significant (approx- 
imately 500%). 

If g(t) denotes the total population of PrPsc, then it evolves eventually as g(t) ~ g(0) exp ( — 
A(V)t) from the ground healthy state [25] in the exponential phase. Here the initial condition 
~ ^inoculation accounts for the inoculation dose. Therefore log g(t) grows linearly in time 
as a first approximation. If we define mathematically the incubation time as a threshold in that 
macroscopic amount of PrPsc [16] , we rederive the logarithmic relationship between incubation 
time and inoculation dose in the linear regime: 

incubation - ~ log ( ^*™*° ) . (25) 
H v ) V ^inoculation / 

Moreover, the relationship (I25p stresses out the importance of measuring the dependence of 
A(V) with respect to the parameters of the model. Following Figure [2] we suggest to evaluate the 
influence of the tightness of the transconformation peak. Closely looking at Figure EJright) we 
indeed notice that the growth rate — A(V) (right limit of the graph) is slightly larger for the model 
with constant rate t{x) = tq (solid line). Several transconformation rates t(x) - being the same 
basal rate To combined with a more and more concentrated gaussian bell are tested - and the 
corresponding growth rate in the exponential phase — A(V) are computed numerically (Figured]). 
Interestingly, the results exhibit a best compromise around a ~ 0.01 (intermediate concentrations 
of the peak). However, this does not correspond to a bimodal distribution for the polymerization 
profile. Thus, according to this model the optimal conditions for PrPsc accumulation in the 
exponential expansion phase are not those which lead to a bimodal size-distribution of PrPsc 
aggregates as observed by Silveira et al. [19J. 



Prion strain mechanisms. It is a challenging problem to investigate the prion strain phe- 
nomenon under the light of model ([2]) which provides informations on the possible microscopic 
distributions of PrPsc polymers. The prion strain phenomenon is supported by the fact that 
prion infected animals may develop several distinct pathologies, whose clinical and neuropatho- 
logical outcomes can be maintained through several passages in rodents models of prion diseases. 
Each prion strain transmitted in the same congenic host exhibits stable incubation period, clin- 
ical expression, lesions profile and infectious properties. These phenotypic features can be used 
to differentiate between different prion strains [7] . Among them, the most commonly used is the 
incubation period. For most prion strains within the same host, incubation period is determined 
by the rate and pattern of PrPsc accumulation. For our purpose however it will be reduced to 
the time for the total amount of PrPsc to reach a certain threshold. We therefore study by mean 
of direct numerical simulations how parameters influence the dynamics of PrPsc accumulation. 
The values used for these numerical tests are mentioned above. The locus and the thightness 
of the transconformation rate are qualitatively compatible with a bimodal distribution of PrPsc 
polymers in the exponential growth phase (see figure [6] and figure [7] B), as well as in the equi- 
librium phase (data not shown). For the experiments depicted in figure [6] and figure [7] varying 
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Figure 5: Influence of the transconformation tightness (top) The transconformation 
rate with several levels of concentration: t(x) = tq + a<p(a(x — m)), where cp is a gaussian 
function and a = 10 ( -~ 3: ' 2:0 )). (middle) The corresponding eigenfunctions U(V;x). (bottom) 
The exponential growth rate — A(V) (solid lind*)Pas a function of a (logarithmic scale, units for 
— A(V) are day -1 ) and the effective transconformation rate (dashed line) r e g- = J t(x)U(V, x) dx 
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Figure 6: Prion replication with variations in the maximum level of the transconformation rate 
t. Transconformation rate is given byr(x) = 0.001 + H*exp(-10* (x-2) 2 ). Three values for If 
have been tested : H = 0.001 (solid line), H = .01 (dashed line) and H = 0.1 (dotted line). (A) 
Size distribution of r (abscissa = PrPsc aggregates size; ordinate = rate r) (B) Time evolution 
of total PrPsc. The arrow represents time t = 96 days (abscissa = Time (in day) ; ordinate = 
rate r (per day) (C) Normalized PrPsc aggregates size distribution at t = 96 days, corresponding 
to the exponential growth of PrPsc. The distributions are normalized by the total number of 
PrPsc aggregates (abscissa = PrPsc aggregates size; ordinate = PrPsc aggregates number). 



parameters acts upon the kinetics of PrPsc accumulation as follows. The higher the transconfor- 
mation or the fragmentation rates are, the faster PrPsc accumulates, and thus, the shorter the 
incubation time seems to be. In particular this implies that a fast strain could be either instable 
(i.e. can be easily broken) or able to fix and transconform PrPc with great efficiency (or both). 
Interestingly, our simulations also predict changes in the microscopic size-distribution of PrPsc 
aggregates, whose profiles are specific of the given parameters. 



6 Conclusion and perspectives 

Our motivation for this work is to better understand the unexpected bimodal size distribution 
of prions, observed by Silveira et al [19]. Based on the nucleated polymerization model of prion 



A B C 




Figure 7: Prion replication with variations in the fragmentation rate of aggregates. Three value 
for have been tested : j3 = 0.0314 (solid line), (3 = 0.0471 (dashed line) and = 0.0628 (dotted 
line). (A), (B) and (C) : Same as Figured] 
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fibril growth describing the dynamics of PrPc monomeres and PrPsc aggregates, we take into 
account aggregate-size dependent parameters, expanding thus a previous study [15] . We have 
shown numerically that bimodal stationary distribution can occur in the case of a non-uniform 
conversion rate from PrPc into PrPsc. It is worth noticing that experimental distribution of PrP 
may not correspond to a steady state of the system since the biological process is dramatically 
stopped with the animal death. Therefore the considerations about the stationary distribution 
might be not biologically relevant. Nevertheless we also observe bimodal distribution generically 
along the temporal dynamics, a new feature in comparison to [16] . For further modifications 
concerning model ([2]) we suggest to take into account saturation in the fragmentation rate for 
large fibrils, as noticed in |19j . 

Theoretical analysis of the model gives necessary criterion on the parameters for a bimodal 
distribution to occur in the stationary phase. Numerical analysis of the eigenvalue problem 
shows in addition that parameters which lead to a bimodal distribution of PrPsc are not the 
most favourable for PrPsc accumulation in the exponential phase, see Figure EJ It suggests that 
infectious ability concentrated on a restricted range of PrPsc aggregate sizes (14 to 28 molecules 
[19] ) does not consist in a replicative advantage for prion at the molecular level. Therefore, this 
restriction could be caused by a physico-chemical constraint, carried out by the host metabolism 
or by the prion aggregate itself. This could eventually be due to the PrPsc aggregate shape, 
since the most infectious intermediate aggregates have been found to be spherical or ellipsoidal 
whereas the others form fibrilar structures [19] . However, understanding the precise nature of 
the constraint would help to identify new therapeutic targets. 

Furthermore, we investigate the implications of varying the fragmentation and transcon- 
formation rates notably in the biological perspectives of the prion strain phenomenon. Prion 
strains have been initially distinguished by incubation periods and lesion profiles in congenic 
mice [U[35]. Nowadays, a large body of literature suggests that differences between prion strains 
lie in the diversity of structures of PrPsc aggregates that can be stably and faithfully propagated 
[6j EH EH EH ES EH HOI HI] • However, it remains poorly understood how these changes in the 
conformation of PrPsc aggregates can account for their physiopathological effects [37] . Among 
the attempted biochemical characterisation of prion strains, a relationship was found between the 
relative stability values of PrPsc aggregates [37] or level of aggregation [38J and incubation times, 
indicating that less stable prions are more infectious, as judged by their shorter incubation times. 
This is presumably because unstable prions fragment more easily, giving rise to smaller aggregates 
of PrPsc that are more infectious than larger ones (represented by an increased transconforma- 
tion rate in the mathematical model). Our model agrees with these observations, since increasing 
the fragmentation rate (3 leads to a faster PrPsc accumulation. More interestingly, our model 
enables to explore quickly the influence of every elementary mechanism involved in prion repli- 
cation. Whereas changing parameters leads to similar effects on PrPsc accumulation kinetics, 
the resulting size-distribution variations are indeed different. The achievement of experimental 
size-distribution of PrPsc aggregates for many prion strains, as obtained by Silveira et al [19] for 
263K prion strain, could therefore allow to better understand the mechanisms involved in prion 
strain phenomenon, by comparing the experimental distribution and the predicted ones. One 
limit to our approach is that prion strains are not only characterized by a precise incubation time 
but also by a specific cellular tropism. It would be of interest to introduce some mathematical 
formalism taking into account cell heterogenity in the brain and the resulting local influences. 

In addition, transmission of prion diseases between different mammalian species is almost 
systematically less efficient than within a single species. This obstruction has been termed 
Species Barriers. Early studies argue that barrier resides in PrP primary structure difference 



18 



between donor and recipient species [42J. However, this issue has been called into question 
notably by BSE strain capability Indeed, different strains propagated in the same host may 
thus have completely different barriers to another species. Consequently, transmission barrier 
appears to depend on prion strain specificities [H2IH]. This is supported by our model. Indeed, 
we learn that the renormalized bimodal distribution of PrPsc aggregates tends to align along 
the eigenvector associated to the dominant eigenvalue A(V) (solid line in Figure [2|left)) [25]. In 
our analysis, we assumed implicitely that the distribution of polymers is initially proportional 
to this eigenvector (which depends upon the the parameters of the model, and thus on the 
strain- host relationship). Otherwise there is some delay for the microscopic distribution to reach 
the good shape which slows down the process of prion accumulation. This might account for 
the strain adaptation mechanism, where primary inoculation is associated with more prolonged 
incubation period than subsequent passages in the same type of hosts |34j . According to these 
observations and the experimental ones, it should be very interesting to study experimentally 
the distribution of aggregates and infectivity between strains which are able or not to cross many 
species barriers. The distributions of PrPsc aggregates could be used to evaluate the influence of 
fragmentation, transconformation rate or infectivity distribution on the strenght of prion strain 
in crossing species barriers. These data could give a new signature of prion strain. Meanwhile the 
elucidation of structural basis of prion strains, strenght barrier species discrimination by these 
prion strain signature could be useful to predict the potential transmissibility of prion strain to 
human, notably for recent atypical scrapie strain, Bovine Amyloidic Spongiform Encephalopathy 
(BASE) strain and Chronic Wasting Disease strain [45] . 

Thus, our model could be used to investigate essential features of prion strains. The next 
step of our work is to approximate faithfully the inverse problem in order to obtain the size 
dependence of the transconformation rate from the distribution of PrPsc. We therefore would be 
able to determine the most infectious compartment for several prion strains. This knowledge is 
a critical step for experimental approaches of prion infectivity investigation, like PMCA. PMCA 
goal is to quickly synthetize in vitro large amounts of PrPsc starting with minute amounts of 
prions. The yield of this technique is very sensitive to experimental procedures. Then size 
distribution of infectivity could help to optimize strain specific PMCA protocols, notably by 
adaptating sonication steps to fall into the most infectious size compartment of PrPsc aggregates. 

Finally, we have analyse the stability of the steady states in a general framework. The diffi- 
culty arising here is that we cannot reduce the study to a set of Ordinary Differential Equations, 
as it is the case in [16} [18]. Under general assumptions on the coefficients, compatible with 
bimodal distributions, the asymptotic stability of the healthy state (i.e. no prion aggregates) is 
established when the PrPc is low. We also prove that this healthy state is unstable when the 
PrPc production rate is high enough. This is in accordance with the asymptotic stability of the 
non-zero steady state for 'constant coefficients' proved in [16J, even though a perfect dichotomy 
between the two results is left open in the general case. Biologically, these results can be inter- 
preted as the propensity of PrPsc aggregates to give rise to prion disease depends on the amount 
of PrPc. In the more general context of understanding why some amyloids are infectious and 
others are not, it therefore could be useful to investigate the amount of amyloidogenic precursors. 
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